﻿///////////////////////////////////////////////////////////////////////////////
// This implements the simple separable convolution
// use only the global memory
///////////////////////////////////////////////////////////////////////////////

#include <assert.h>
#include "../common/common.h"

#define IMAGE_WIDTH 2048
#define IMAGE_HEIGHT 2048

#define KERNEL_RADIUS 8
#define KERNEL_LENGTH (2 * KERNEL_RADIUS + 1)

#define ROW_BLOCKDIM_X 16
#define ROW_BLOCKDIM_Y 16
#define COLUMN_BLOCKDIM_X 16
#define COLUMN_BLOCKDIM_Y 16

/*==============================================================================*/
// CPU Separable Convolution
/*==============================================================================*/

void convolutionColumnCPU(float *h_Dst, float *h_Src, float *h_Kernel,
						  int imageW, int imageH, int kernelR)
{
	for(int y = 0; y < imageH; y++)
		for(int x = 0; x < imageW; x++){
			float sum = 0;
			for(int k = -kernelR; k <= kernelR; k++){
				int dy = y + k;
				if(dy >= 0 && dy < imageH)
					sum += h_Src[dy * imageW + x] * h_Kernel[kernelR - k];
			}
			h_Dst[y * imageW + x] = sum;
		}
}

void convolutionRowCPU(float *h_Dst, float *h_Src, float *h_Kernel,
					   int imageW, int imageH, int kernelR)
{
	for(int y = 0; y < imageH; y++)
		for(int x = 0; x < imageW; x++){
			float sum = 0;
			for(int k = -kernelR; k <= kernelR; k++){
				int dx = x + k;
				if(dx >= 0 && dx < imageW)
					sum += h_Src[y * imageW + dx] * h_Kernel[kernelR - k];
			}
			h_Dst[y * imageW + x] = sum;
		}
}


/*==============================================================================*/
// GPU Separable Convolution
/*==============================================================================*/

__global__ void convolutionColumnGPUKernel( float *d_Dst, float *d_Src, float *c_VerticalKernel)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	float sum = 0;

	for(int k = -KERNEL_RADIUS; k <= KERNEL_RADIUS; k++){
		int dy = y + k;
		if(dy >= 0 && dy < IMAGE_HEIGHT)
			sum += d_Src[dy * IMAGE_WIDTH + x] * c_VerticalKernel[KERNEL_RADIUS - k];
	}
	d_Dst[y * IMAGE_WIDTH + x] = sum;
}

__global__ void convolutionRowGPUKernel( float *d_Dst, float *d_Src, float *c_HorizontalKernel)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	float sum = 0;

	for(int k = -KERNEL_RADIUS; k <= KERNEL_RADIUS; k++){
		int dx = x + k;
		if(dx >= 0 && dx < IMAGE_WIDTH)
			sum += d_Src[y * IMAGE_WIDTH + dx] * c_HorizontalKernel[KERNEL_RADIUS - k];
	}
	d_Dst[y * IMAGE_WIDTH + x] = sum;
}

/*==============================================================================*/
// main() function
/*==============================================================================*/

int main()
{
	float	*h_VerticalKernel,
			*h_HorizontalKernel,
			*h_Input,
			*h_Buffer,
			*h_OutputCPU,
			*h_OutputGPU;
	float	*d_Input,
			*d_Buffer,
			*d_Output;
	float	*c_VerticalKernel,
			*c_HorizontalKernel;
	float	elapsedTime;
	cudaEvent_t start, stop;
	
	// Cấp phát bộ nhớ trên CPU
	h_VerticalKernel = (float *)malloc(KERNEL_LENGTH * 
		sizeof(float));
	h_HorizontalKernel = (float *)malloc(KERNEL_LENGTH *
		sizeof(float));
	h_Input = (float *)malloc(IMAGE_WIDTH * IMAGE_HEIGHT *
		sizeof(float));
	h_Buffer = (float *)malloc(IMAGE_WIDTH * IMAGE_HEIGHT * 
		sizeof(float));
	h_OutputCPU = (float *)malloc(IMAGE_WIDTH * IMAGE_HEIGHT *
		sizeof(float));
	h_OutputGPU = (float *)malloc(IMAGE_WIDTH * IMAGE_HEIGHT *
		sizeof(float));

	// Phát sinh dữ liệu ngẫu nhiên
	// cho mảng ban đầu và 2 vector bộ lọc
	srand(200);
	for(unsigned int i = 0; i < KERNEL_LENGTH; i++){
		h_VerticalKernel[i] = (float)(rand() % 16);
		h_HorizontalKernel[i] = (float)(rand() % 16);
	}
	for(unsigned int i = 0; i < IMAGE_WIDTH * IMAGE_HEIGHT; i++)
		h_Input[i] = (float)(rand() % 16);

	// Ghi nhận thời điểm bắt đầu của đoạn mã GPU
	SAFE_CALL( cudaEventCreate( &start ) );
	SAFE_CALL( cudaEventCreate( &stop ) );
	SAFE_CALL( cudaEventRecord( start, 0 ) );

	// Cấp phát bộ nhớ GPU
	SAFE_CALL( cudaMalloc(&d_Input, IMAGE_WIDTH * IMAGE_HEIGHT * sizeof(float)) );
	SAFE_CALL( cudaMalloc(&d_Buffer, IMAGE_WIDTH * IMAGE_HEIGHT * sizeof(float)) );
	SAFE_CALL( cudaMalloc(&d_Output, IMAGE_WIDTH * IMAGE_HEIGHT * sizeof(float)) );

	// Cấp phát bộ nhớ cho 2 vector bộ lọc
	SAFE_CALL( cudaMalloc(&c_HorizontalKernel, KERNEL_LENGTH * sizeof(float)) );
	SAFE_CALL( cudaMalloc(&c_VerticalKernel, KERNEL_LENGTH * sizeof(float)) );

	// chép 2 vector bộ lọc từ CPU sang GPU
	SAFE_CALL ( cudaMemcpy(c_HorizontalKernel, h_HorizontalKernel, 
		KERNEL_LENGTH * sizeof(float), cudaMemcpyHostToDevice) );
	SAFE_CALL ( cudaMemcpy(c_VerticalKernel, h_VerticalKernel, 
		KERNEL_LENGTH * sizeof(float), cudaMemcpyHostToDevice) );

	// chép dữ liệu input từ CPU sang GPU
	SAFE_CALL ( cudaMemcpy(d_Input, h_Input, 
		IMAGE_WIDTH * IMAGE_HEIGHT * sizeof(float), cudaMemcpyHostToDevice) );

	// Chuẩn bị tham số cấu hình cho các hàm kernel
	assert(IMAGE_WIDTH % ROW_BLOCKDIM_X == 0);
	assert(IMAGE_HEIGHT % ROW_BLOCKDIM_Y == 0);
	assert(IMAGE_WIDTH % COLUMN_BLOCKDIM_X == 0);
	assert(IMAGE_HEIGHT % COLUMN_BLOCKDIM_Y == 0);

	// Gọi bước đầu tiên của thuật toán trên GPU
	dim3 threads(COLUMN_BLOCKDIM_X, COLUMN_BLOCKDIM_Y);
	dim3 blocks(IMAGE_WIDTH / COLUMN_BLOCKDIM_X, IMAGE_HEIGHT / COLUMN_BLOCKDIM_Y);
	convolutionColumnGPUKernel<<<blocks, threads>>>(d_Buffer, d_Input, c_VerticalKernel);

	// Bước thứ 2 của thuật toán
	threads = dim3(ROW_BLOCKDIM_X, ROW_BLOCKDIM_Y);
	blocks = dim3(IMAGE_WIDTH / ROW_BLOCKDIM_X, IMAGE_HEIGHT / ROW_BLOCKDIM_Y);
	convolutionRowGPUKernel<<<blocks, threads>>>(d_Output, d_Buffer, c_HorizontalKernel);

	// Chép kết quả từ GPU sang CPU
	SAFE_CALL (cudaMemcpy(h_OutputGPU, d_Output, 
		IMAGE_WIDTH * IMAGE_HEIGHT * sizeof(float), cudaMemcpyDeviceToHost) );

	// Ghi nhận thời điểm kết thúc
	// và hiển thị thời gian thực thi
	SAFE_CALL( cudaEventRecord( stop, 0 ) );
	SAFE_CALL( cudaEventSynchronize( stop ) );
	SAFE_CALL( cudaEventElapsedTime( &elapsedTime, start, stop ) );

	printf( "Thoi gian thuc thi %3.1f ms\n", elapsedTime );
	SAFE_CALL( cudaEventDestroy( start ) );
	SAFE_CALL( cudaEventDestroy( stop ) );

	// Gọi lần lượt 2 hàm tính tích chập trên CPU
	convolutionColumnCPU(h_Buffer, h_Input, h_VerticalKernel,
		IMAGE_WIDTH, IMAGE_HEIGHT, KERNEL_RADIUS);
	convolutionRowCPU(h_OutputCPU, h_Buffer, h_HorizontalKernel,
		IMAGE_WIDTH, IMAGE_HEIGHT, KERNEL_RADIUS);

	// Kiểm tra kết quả giữa h_OutputCPU và h_OutputGPU
	double sum = 0, delta = 0;
	for(unsigned int i = 0; i < IMAGE_WIDTH * IMAGE_HEIGHT; i++){
		delta += (h_OutputGPU[i] - h_OutputCPU[i]) * (h_OutputGPU[i] - h_OutputCPU[i]);
		sum   += h_OutputCPU[i] * h_OutputCPU[i];
	}
	double L2norm = sqrt(delta / sum);
	printf("Relative L2 norm: %E\n\n", L2norm);
	printf((L2norm < 1e-6) ? "OK\n\n" : "FAILED\n\n");

	// giải phóng bộ nhớ
	free(h_OutputCPU);
	free(h_Buffer);
	free(h_Input);
	free(h_VerticalKernel);
	free(h_HorizontalKernel);
	free(h_OutputGPU);

	SAFE_CALL( cudaFree(d_Input) );
	SAFE_CALL( cudaFree(d_Buffer) );
	SAFE_CALL( cudaFree(d_Output) );
	SAFE_CALL( cudaFree(c_HorizontalKernel) );
	SAFE_CALL( cudaFree(c_VerticalKernel) );

	return 0;
}
